function [Z,Tris,Points] = fe_solution_bb(V, T, x, FE_Order, triangulars, d_plot)
%function [Z,Tris,Points] = fe_solution_bb(V,T,x,FE_Order,triangulars,d_plot)
% this is the version for single degrees.

m = (d_plot+1)*(d_plot+2)/2;
n_dof_per_elem = (FE_Order + 1)*(FE_Order + 2)/2;
[I,J,K] = indices(d_plot);
cell_per_tri = size(triangulars,1);
mat = vdm22_bb(FE_Order, d_plot);

n_elem = size(T,1);
Tris = zeros(cell_per_tri*n_elem,3);
Points = zeros(m*n_elem,2);
Z = zeros(m*n_elem,1);
pos = 0;
for tri = 1:n_elem
    Points((tri-1)*m+1:tri*m,:) = (I*V(T(tri,1),:)+J*V(T(tri,2),:)+K*V(T(tri,3),:))/d_plot;
    Z((tri-1)*m+1:tri*m) = mat * x(:,tri);
    Tris((tri-1)*cell_per_tri+1:tri*cell_per_tri,:) = triangulars + (tri-1)*m;
    
    pos = pos + n_dof_per_elem;
end

end

function Mat= vdm22_bb(d_col,d_row)
% d_row represent the collotation points degree;
% d_col represent the degree of bform basis
m_row = (d_row + 1)*(d_row + 2)/2;
m_col = (d_col + 1)*(d_col + 2)/2;
[I_row,J_row,K_row] = indices(d_row);
[I_col,J_col,K_col] = indices(d_col);
IM_row = diag(I_row)*ones(m_row,m_col);
JM_row = diag(J_row)*ones(m_row,m_col);
KM_row = diag(K_row)*ones(m_row,m_col);
IM_col = diag(I_col)*ones(m_col,m_row);
JM_col = diag(J_col)*ones(m_col,m_row);
KM_col = diag(K_col)*ones(m_col,m_row);
Mat = (IM_row/d_row).^(IM_col').*(JM_row/d_row).^(JM_col').*(KM_row/d_row).^(KM_col');
IF = gamma(I_col+1);
JF = gamma(J_col+1);
KF = gamma(K_col+1);
A = factorial(d_col)*ones(m_row,m_col)*diag(1./(IF.*JF.*KF));
Mat = A.*Mat;
end